/**
 Application.cxx
 \file Definition of class Application
 
 \author Thomas Burton
 \date 12/6/11
 \copyright 2011 BNL
 */

#include "tppmc/Application.h"

#include <cmath>

#include <TROOT.h>
#include <TCollection.h>
#include <TNamed.h>

#include <TClonesArray.h>
#include <TDirectory.h>
#include <TMath.h>
#include <TPythia6.h>
#include <TTree.h>

#include "tppmc/DistributionQaFigures.h"
#include "tppmc/Distributions.h"
#include "tppmc/Enumerations.h"
#include "tppmc/Event.h"
#include "tppmc/FragmentationFunction.h"
#include "tppmc/PartonDistribution.h"

namespace tppmc {

   Application::Application() {
      Splash(std::cout);
   }

   /** \todo add try...catch in case of bad configuration */
   Application::Application(Configuration& config)
   : mConfig(config)
   , mTree(NULL)
   , mEvent(NULL) {
      Splash(std::cout);
      Initialise(config);
   }

   Application::~Application() {
      TIter next(gROOT->GetListOfFunctions());
      // Most pointers are held by smart pointers.
      // TTree is owned by the output file and is deleted when
      // the file closes, so we don't delete it.
      if(mEvent) {
         delete mEvent;
         mEvent = NULL;
      } // if
   }

   void Application::Splash(std::ostream& os) const {
      os << std::endl <<
      "                 TPPMC" << std::endl <<
      "Transverse Polarized Proton Monte Carlo" << std::endl <<
      "              Version " << smVersionMajor << "." << smVersionMinor <<
      std::endl << std::endl <<
      "Monte Carlo Generator for asymmetry studies\n" <<
      "   in polarized proton-proton collisions\n\n" <<
      "      http://code.google.com/p/tppmc/" <<
      std::endl;
   }

   bool Application::Initialise(Configuration& config) {
      if(not config.Verify()) {
         return false;
      } // if
      // Suppress printing of messages about pages added to postscript,
      // new canvases etc
      gErrorIgnoreLevel = 5000;
      
      mConfig = config;
      std::cout << "\nHere is your run configuration:\n" << std::endl;
      config.Speak(std::cout);
      
      mRandom.reset(new TRandom3(0));
      
      mPdfs.reset(new Distributions(config));
      std::cout << "The distribution functions to be used are:" << std::endl;
      mPdfs->Speak();
      
      mTreeFile.reset(new TFile(config.outfile.c_str(), "RECREATE"));
      if(mTree) {
         delete mTree;
      } // if
      mTree = new TTree("events", "TPPMC event tree");
      if(mEvent) {
         delete mEvent;
      } // if
      mEvent = new Event(0., config.EtaMin(), config.EtaMax(),
                         config.MinHadronPt());
      mTree->Branch("event", &mEvent, 32000, 99);
      
      if(config.GenerateQaFigures()) {
         std::cout << "Generating qa figures to " << config.QaFiguresFileName() << std::endl;
         TDirectory* dir = gDirectory; // Remember initial directory
         
         mQaFile.reset(new TFile(config.QaFiguresFileName().c_str(),
                                 "RECREATE"));
         
         tppmc::DistributionQaFigures figures(NULL);
         figures.Fill(*mPdfs);
         figures.Write();
         
         // Restore the initial directory as the current directory
         gDirectory = dir;
         if(dir) dir->cd();
         std::cout << "Finished QA figures" << std::endl;
      } // if
      // Creating a TStopwatch also starts it, so we do this last -
      // we immediately go into the event loop after Initialise() finishes.
      mTimer.reset(new TStopwatch);
      return true;
   }

   bool Application::ProcessEvent() {
      /**
       \todo Protect status message from printing multiple times for a given event
       count if there are cuts on events (eta, pt).
       */
      if(mTree->GetEntries() % 100 == 0 and mTree->GetEntries() not_eq 0) {
         // Print the time since the timer began
         std::cout << "Generated " << mTree->GetEntries() << " events in "
         << mTimer->RealTime() // Stops the timer
         << " seconds" << std::endl;
         // Restart the timer
         mTimer->Start(true);
      } // if
      // Generate a vertex for this event 
      if((mConfig.machine_type not_eq FNAL_E704_PROTON) and
         (mConfig.machine_type not_eq FNAL_E704_PBAR)) {
         double zvtx(0.);
         do {
            zvtx = mRandom->Gaus(0.0, 23.);
         } while(std::fabs(zvtx) > 30.);
      } // if
      // Generate a pythia event and get the particles
      TPythia6::Instance()->GenerateEvent();
      // Rotate the event by a phi angle
      // (pythia generates scattered partons at fixed phi)
      mEvent->SetRotation(mRandom->Uniform(-TMath::Pi(), TMath::Pi()));
      mEvent->Update(*mPdfs);
      // Fill the tree if we have hadrons.
      // If we are not fragmentating (function = 0) always fill, because
      // we never have hadrons but we want to save the parton info.
      if(mEvent->Hadrons().size() > 0 or mConfig.FragmentationFunction() == 0) {
         mTree->Fill();
      } // if
      return true;
   }

   bool Application::Close() {
       mPdfs->Unpolarised()->Write("unpol");
       mPdfs->Transversity()->Write("transversity");
       mPdfs->Sivers()->Write("sivers");
       mPdfs->Fragmentation()->Write("frag");
       return mTreeFile->Write() > 0;
   }

   bool Application::EventLoop(int n) {
      if(n < 1) {
         n = mConfig.NEvents();
      } // if
      while(mTree->GetEntries() < n) {
         ProcessEvent();
      } // while
      return mTree->GetEntries() == n;
   }
} // namespace tppmc
